##########################################
### Analysis: MSA Models #################
##########################################

## Purpose: This script includes the models for the
## MSA analysis.
## Data In: 
## 1) data_files/datamsa_4_30_2022.rds
## Data Out: 
## 1) plots/AJS/local_year_alternative
## 2) plots/AJS/segregation_year
## 3) plots/AJS/local_party_time.pdf
## 4) plots/AJS/local_ses_time.pdf
## Software Versions:

## R version 4.1.2

library(RColorBrewer) ## 1.1.2
library(tidyverse) ## 1.3.1
library(viridis) ## 0.6.2
library(MASS) ## 7.3.54

## set working directory
## to replication package

data.MSA <- readRDS("data_files/datamsa_5_6_2025.rds")


## colorblind palette for figures
cbbPalette <- c("#000000", "#E69F00", 
                "#56B4E9", "#009E73", 
                "#F0E442", "#0072B2", 
                "#D55E00", "#CC79A7")


'%ni%' <- Negate("%in%")

## this function calculates s.e. and point estimates
## from parametric bootstrap draws used below

estimates.fct <- function(frame){
  est <- as.data.frame(matrix(NA, ncol = 4, nrow = ncol(frame)))
  colnames(est) <- c("Est", "Lower", "Upper", "Model")
  est$Est <- base::apply(frame, 2, function(x){mean(x)})
  est$Lower <- base::apply(frame, 2, function(x){quantile(x, .025)})
  est$Upper <- base::apply(frame, 2, function(x){quantile(x, .975)})
  est$Model <- colnames(frame)  
  return(est)
}

#############################################
### Creating Year Buckets ###################
############################################


## creating year buckets
data.MSA$year_bucket <- ifelse(data.MSA$year <= 1988,
                               "1987, 1988",
                               ifelse(data.MSA$year <= 1999,
                                      "1992, 1997, 1999",
                                      ifelse(data.MSA$year <= 2007,
                                             "2003, 2007",
                                             "2009, 2012")))


################################
#### Year:Local Interactions ###
################################

## In this section we investigate
## how inequality perceptions vary by
## differing levels of our MSA
## local inequality variables. We include
## plots from this invesitgation in the
## main text, Figure 4 and in the SI,
## section E. 


msa1 <- lm(inequality.num ~ factor(year_bucket)*gini.std + party.harm + race2 + ses + cohort2 + gender,
           data = data.MSA)
msa2 <- lm(inequality.num ~ factor(year_bucket)*segregation.std + party.harm + race2 + ses + cohort2 + gender,
           data = data.MSA)
msa3 <- lm(inequality.num ~ factor(year_bucket)*incexpose.std + party.harm + race2 + ses + cohort2 + gender,
           data = data.MSA)
msa4 <- lm(inequality.num ~ factor(year_bucket)*iGINIhouseVALmed.std + party.harm + race2 + ses + cohort2 + gender,
           data = data.MSA)
msa5 <- lm(inequality.num ~ factor(year_bucket)*iGINIhouseVALq4.std + party.harm + race2 + ses + cohort2 + gender,
           data = data.MSA)
msa6 <- lm(inequality.num ~ factor(year_bucket)*povertyrate.std + party.harm + race2 + ses + cohort2 + gender,
           data = data.MSA)
msa7 <- lm(inequality.num ~ factor(year_bucket)*rmd.std + party.harm + race2 + ses + cohort2 + gender,
           data = data.MSA)

msas <- list(msa1,
             msa2,
             msa3,
             msa4,
             msa5,
             msa6,
             msa7)


### creating small multiples plot over time
### Note: this is the same code
### introduced in 03_analysis_full_data_ev.R
### see that code file for an overview
### of the code structure 
 
estimates_msas <- list()

for(x in 1:length(msas)){
  
  vcov <- vcov(msas[[x]])
  betas <- msas[[x]]$coefficients
  
  set.seed(08540)
  beta.draws <- MASS::mvrnorm(10000, mu = betas, Sigma = vcov)
  
  
  ## create matrix number of estimates by number of covariates
  ## number of estimates: 3 (we're doing estimates for
  ## mean, 1 standard deviation above and below mean for each
  ## estimates) x number of year (4) = 12
  
  matrix <- matrix(0, nrow = 12, 
                   ncol = length(beta.draws[1,]))
  matrix[,1] <- 1
  colnames(matrix) <- names(beta.draws[1,])
  rownames(matrix) <- 1:nrow(matrix)
  
  ## year - different year for each level
  ## main year effects:
  
  ## for each level:
  ## each level gets a row of zeros (for first 1966-1970)
  ## and then a single year 
  for(i in 2:4){
    matrix[i, i] <- 1
    matrix[i+4, i] <- 1
    matrix[i+8, i] <- 1
  } 
  
  
  ## geo variables
  matrix[,5] <- c(rep(-1, 4),
                  rep(0, 4),
                  rep(1, 4))
  
  ## other variables set constant:
  matrix[,"race2Yes"] <- 1
  matrix[,"ses3"] <- 1
  matrix[,"gendermale"] <- 1
  ## Boomers
  
  ## geo-year interaction
  ## first 5: -1 
  ## second 5: 0
  ## third 5: 1
  for(i in 2:4){
    matrix[i, i + 17] <- -1
    matrix[i+4, i + 17] <- 0
    matrix[i+8, i + 17] <- 1
  }
  

  sim <- base::apply(beta.draws, 1, function(x){
    vector <- matrix %*% x
    return(vector)
  })
  sim <- t(sim)
  
  estimates_msas[[x]] <- estimates.fct(sim)
  
  print(x)
}

## cleaning estimate files
for(i in 1:length(estimates_msas)){
  estimates_msas[[i]]$year <- levels(factor(data.MSA$year_bucket))
  estimates_msas[[i]]$year <- dplyr::recode(estimates_msas[[i]]$year,
                                            `1987, 1988` = "1987",
                                            `1992, 1997, 1999` = "1995",
                                            `2003, 2007` = "2005",
                                            `2009, 2012` = "2010")
  estimates_msas[[i]]$year <- as.Date(paste0(estimates_msas[[i]]$year,
                                             "-01-01"))
  estimates_msas[[i]]$std_deviations <- c(rep(-1, 4),
                                          rep(0, 4),
                                          rep(1, 4))
  estimates_msas[[i]]$variable <- names(msas[[i]]$coefficients)[5]
}

estimates_msas <- bind_rows(estimates_msas)



estimates_msas$variable <- dplyr::recode(
  estimates_msas$variable,
  `gini.std` = "Local Gini",
  `segregation.std` = "Segregation",
  `incexpose.std` = "Income Exposure",
  `iGINIhouseVALmed.std` = "Wealth Gini: Median",
  `iGINIhouseVALq4.std` = "Wealth Gini: Q4",
  `povertyrate.std` = "Poverty Rate",
  `rmd.std` = "Relative Mean Deprivation")


ggplot(estimates_msas[estimates_msas$std_deviations != 0 &
                        estimates_msas$variable != "Income Exposure" &
                        estimates_msas$variable != "Segregation", ],
       aes(x = year,
           color = factor(std_deviations),
           y = Est,
           ymin = Lower,
           ymax = Upper)) +
  facet_wrap(~ variable,
             ncol = 2) +
  geom_point(position = position_dodge(2 * 365)) +
  geom_errorbar(width = 0,
                position = position_dodge(2 * 365)) +
  geom_line() +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(x = NULL,
       y = "Estimated Percent of Respondents\nWho Think 'Rich Getting Richer'",
       color = "Standard Deviations from\nSurvey Mean") +
  scale_color_manual(values = cbbPalette[c(2,3)])
## figures/local_year_alternative
## Figure E.1 in the SI
## 8x7

## for main section of paper
## Figure 4
estimates_msas$variable <- dplyr::recode(estimates_msas$variable,
                                         `Segregation` = "Rank Order\nInformation Theory Index")

ggplot(estimates_msas[estimates_msas$std_deviations != 0 &
                        (estimates_msas$variable == "Rank Order\nInformation Theory Index" |
                           estimates_msas$variable == "Income Exposure"), ],
       aes(x = year,
           color = factor(std_deviations),
           lty = factor(std_deviations),
           shape = factor(std_deviations),
           y = Est,
           ymin = Lower,
           ymax = Upper)) +
  facet_wrap(~ variable) +
  geom_line() +
  geom_point(position = position_dodge(2 * 365)) +
  geom_errorbar(width = 0,
                position = position_dodge(2 * 365)) +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(0,1)) +
  labs(x = NULL,
       y = "Estimated Percent of Respondents\nWho Think 'Rich Getting Richer'",
       color = "Standard Deviations\nfrom Survey Mean",
       lty = "Standard Deviations\nfrom Survey Mean",
       shape = "Standard Deviations\nfrom Survey Mean") +
  scale_color_manual(values = cbbPalette[c(2,3)])
## figures/segregation_year
## 5x7
## bw version: figures/segregation_year_bw

########################################
#### local:year party interactions ###########
########################################

## In this section we investigate
## how inequality perceptions vary
## by levels of local MSA inequality
## and partisanship. These results
## are included in the SI, section E. 


msa1 <- lm(inequality.num ~ factor(year_bucket)*incexpose.std*party.harm + cohort2 + ses + race2 + gender,
           data = data.MSA)
msa2 <- lm(inequality.num ~ factor(year_bucket)*segregation.std*party.harm  + cohort2 + ses +race2 + gender,
           data = data.MSA)
msa3  <- lm(inequality.num ~ factor(year_bucket)*segregation.std*ses  + cohort2 + party.harm +race2 + gender,
            data = data.MSA)
msa4  <- lm(inequality.num ~ factor(year_bucket)*incexpose.std*ses  + cohort2 + party.harm +race2 + gender,
            data = data.MSA)

msas <- list(msa1,
             msa2,
             msa3,
             msa4)

### creating small multiples plot over time
estimates_msas <- list()

for(x in 1:length(msas)){
  
  vcov <- vcov(msas[[x]])
  betas <- msas[[x]]$coefficients
  
  set.seed(08540)
  beta.draws <- MASS::mvrnorm(10000, mu = betas, Sigma = vcov)
  
  ## create matrix number of estimates by number of covariates
  ## number of estimates: 3 (we're doing estimates for
  ## mean, 1 standard deviation above and below mean for each
  ## estimates) x number of year (4) x party (3) = 36
  
  ## for ses (3 and 4)
  ## it's 3 x 4 x 5 = 60
  
  if(x %in% c(1,2)){
    matrix <- matrix(0, nrow = 36, 
                     ncol = length(beta.draws[1,]))
    matrix[,1] <- 1
    colnames(matrix) <- names(beta.draws[1,])
    rownames(matrix) <- 1:nrow(matrix)
    
    ## year - different year for each level
    ## main year effects:
    
    ## for each level:
    ## each level gets a row of zeros (for first 1966-1970)
    ## and then a single year 
    ## Democrat - 3 sets of 4 years 
    ## Republican - 3 sets of 4 years
    ## Indep/other - 3 sets of 4 years
    for(i in 2:4){
      matrix[i, i] <- 1
      matrix[i+4, i] <- 1
      matrix[i+8, i] <- 1
      matrix[i+12, i] <- 1
      matrix[i+16, i] <- 1
      matrix[i+20, i] <- 1
      matrix[i+24, i] <- 1
      matrix[i+28, i] <- 1
      matrix[i+32, i] <- 1
    } 
    
    ## geo variables
    matrix[,5] <- rep(c(rep(-1, 4),
                        rep(0, 4),
                        rep(1, 4)), 3)
    
    matrix[,"party.harmDemocrat"] <- c(rep(1,12),
                                       rep(0,24))
    matrix[,"party.harmIndep/Other"] <- c(rep(0,12),
                                          rep(0,12),
                                          rep(1,12))
    ## boomers
    matrix[,"ses3"] <- 1
    matrix[,"race2Yes"] <- 1
    matrix[,"gendermale"] <- 1
    
    ## geo-year interaction
    ## first 4: -1 
    ## second 5: 0
    ## third 5: 1
    for(i in 2:4){
      matrix[i, i + 17] <- -1
      matrix[i+4, i + 17] <- 0
      matrix[i+8, i + 17] <- 1
      matrix[i+12, i + 17] <- -1
      matrix[i+16, i + 17] <- 0
      matrix[i+20, i + 17] <- 1
      matrix[i+24, i + 17] <- -1
      matrix[i+28, i + 17] <- 0
      matrix[i+32, i + 17] <- 1
    }
    
    ## year:party interaction
    for(i in 2:4){
      matrix[i, i + 20] <- 1
      matrix[i +4, i + 20] <- 1
      matrix[i + 8, i + 20] <- 1
      matrix[i + 24, i + 23] <- 1
      matrix[i + 28, i + 23] <- 1
      matrix[i + 32, i + 23] <- 1
    }
    
    ## gini:party
    if(x == 1){
      matrix[1:12, "incexpose.std:party.harmDemocrat"] <- c(rep(-1, 4),
                                                       rep(0, 4),
                                                       rep(1, 4))
      matrix[25:36, "incexpose.std:party.harmIndep/Other"] <- c(rep(-1, 4),
                                                           rep(0, 4),
                                                           rep(1, 4))
      ## year: gini:party interaction
      matrix[c(2,6,10), "factor(year_bucket)1992, 1997, 1999:incexpose.std:party.harmDemocrat"] <- c(-1,0,1)
      matrix[c(3,7,11), "factor(year_bucket)2003, 2007:incexpose.std:party.harmDemocrat"] <- c(-1,0,1)
      matrix[c(4,8,12), "factor(year_bucket)2009, 2012:incexpose.std:party.harmDemocrat"] <- c(-1,0,1)
      matrix[c(26,30,34), "factor(year_bucket)1992, 1997, 1999:incexpose.std:party.harmIndep/Other"] <- c(-1,0,1)
      matrix[c(27,31,35), "factor(year_bucket)2003, 2007:incexpose.std:party.harmIndep/Other"] <- c(-1,0,1)
      matrix[c(28,32,36), "factor(year_bucket)2009, 2012:incexpose.std:party.harmIndep/Other"] <- c(-1,0,1)
      
    } else {
      matrix[1:12, "segregation.std:party.harmDemocrat"] <- c(rep(-1, 4),
                                                              rep(0, 4),
                                                              rep(1, 4))
      matrix[25:36, "segregation.std:party.harmIndep/Other"] <- c(rep(-1, 4),
                                                                  rep(0, 4),
                                                                  rep(1, 4))
      ## year: gini:party interaction
      matrix[c(2,6,10), "factor(year_bucket)1992, 1997, 1999:segregation.std:party.harmDemocrat"] <- c(-1,0,1)
      matrix[c(3,7,11), "factor(year_bucket)2003, 2007:segregation.std:party.harmDemocrat"] <- c(-1,0,1)
      matrix[c(4,8,12), "factor(year_bucket)2009, 2012:segregation.std:party.harmDemocrat"] <- c(-1,0,1)
      matrix[c(26,30,34), "factor(year_bucket)1992, 1997, 1999:segregation.std:party.harmIndep/Other"] <- c(-1,0,1)
      matrix[c(27,31,35), "factor(year_bucket)2003, 2007:segregation.std:party.harmIndep/Other"] <- c(-1,0,1)
      matrix[c(28,32,36), "factor(year_bucket)2009, 2012:segregation.std:party.harmIndep/Other"] <- c(-1,0,1)
      
    }
    
    
  } else {
    
    matrix <- matrix(0, nrow = 60, 
                     ncol = length(beta.draws[1,]))
    matrix[,1] <- 1
    colnames(matrix) <- names(beta.draws[1,])
    rownames(matrix) <- 1:nrow(matrix)
    
    ## year - different year for each level
    ## main year effects:
    
    ## for each level:
    ## each level gets a row of zeros (for first 1966-1970)
    ## and then a single year 
    ## ses1 - 3 sets of 4 years 
    ## ses2 - 3 sets of 4 years
    ## ses3 - 3 sets of 4 years
    ## ses4 - 3 sets of 4 years
    ## ses5 - 3 sets of 4 years
    for(i in 2:4){
      matrix[i, i] <- 1
      matrix[i+4, i] <- 1
      matrix[i+8, i] <- 1
      matrix[i+12, i] <- 1
      matrix[i+16, i] <- 1
      matrix[i+20, i] <- 1
      matrix[i+24, i] <- 1
      matrix[i+28, i] <- 1
      matrix[i+32, i] <- 1
      matrix[i+36, i] <- 1
      matrix[i+40, i] <- 1
      matrix[i+44, i] <- 1
      matrix[i+48, i] <- 1
      matrix[i+52, i] <- 1
      matrix[i+56, i] <- 1
    } 
    
    ## geo variables
    matrix[,5] <- rep(c(rep(-1, 4),
                        rep(0, 4),
                        rep(1, 4)), 5)
    
    matrix[,"party.harmDemocrat"] <- 0
    matrix[,"party.harmIndep/Other"] <- 0
    matrix[,"race2Yes"] <- 1
    matrix[, "gendermale"] <- 1
    
    ## boomers
    matrix[,"ses2"] <- c(rep(0, 12),
                         rep(1, 12),
                         rep(0, 36))
    matrix[,"ses3"] <-  c(rep(0, 24),
                          rep(1, 12),
                          rep(0, 24))
    matrix[,"ses4"] <-  c(rep(0, 36),
                          rep(1, 12),
                          rep(0, 12))
    matrix[,"ses5"] <-  c(rep(0, 48),
                          rep(1, 12))
    
    
    ## geo-year interaction
    ## first 4: -1 
    ## second 4: 0
    ## third 4: 1
    for(i in 2:4){
      matrix[i, i + 17] <- -1
      matrix[i+4, i + 17] <- 0
      matrix[i+8, i + 17] <- 1
      matrix[i+12, i + 17] <- -1
      matrix[i+16, i + 17] <- 0
      matrix[i+20, i + 17] <- 1
      matrix[i+24, i + 17] <- -1
      matrix[i+28, i + 17] <- 0
      matrix[i+32, i + 17] <- 1
      matrix[i+36, i + 17] <- - 1
      matrix[i+40, i + 17] <- 0
      matrix[i+44, i + 17] <- 1
      matrix[i+48, i + 17] <- -1
      matrix[i+52, i + 17] <- 0
      matrix[i+56, i + 17] <- 1
    }
    
    ## year:ses interaction
    for(i in 2:4){
      matrix[i + 12, i + 20] <- 1
      matrix[i + 16, i + 20] <- 1
      matrix[i + 20, i + 20] <- 1
      matrix[i + 24, i + 23] <- 1
      matrix[i + 28, i + 23] <- 1
      matrix[i + 32, i + 23] <- 1
      matrix[i + 36, i + 26] <- 1
      matrix[i + 40, i + 26] <- 1
      matrix[i + 44, i + 26] <- 1
      matrix[i + 48, i + 29] <- 1
      matrix[i + 52, i + 29] <- 1
      matrix[i + 56, i + 29] <- 1
    }
    
    ## gini:ses
    if(x == 4){
      matrix[13:24, "incexpose.std:ses2"] <- c(rep(-1, 4),
                                          rep(0, 4),
                                          rep(1, 4))
      matrix[25:36, "incexpose.std:ses3"] <- c(rep(-1, 4),
                                          rep(0, 4),
                                          rep(1, 4))
      matrix[37:48, "incexpose.std:ses4"] <- c(rep(-1, 4),
                                          rep(0, 4),
                                          rep(1, 4))
      matrix[49:60, "incexpose.std:ses5"] <- c(rep(-1, 4),
                                          rep(0, 4),
                                          rep(1, 4))
      
      ## year: gini:ses interaction
      matrix[c(14,18,22), "factor(year_bucket)1992, 1997, 1999:incexpose.std:ses2"] <- c(-1,0,1)
      matrix[c(15,19,23), "factor(year_bucket)2003, 2007:incexpose.std:ses2"] <- c(-1,0,1)
      matrix[c(16,20,24), "factor(year_bucket)2009, 2012:incexpose.std:ses2"] <- c(-1,0,1)
      matrix[c(26,30,34), "factor(year_bucket)1992, 1997, 1999:incexpose.std:ses3"] <- c(-1,0,1)
      matrix[c(27,31,35), "factor(year_bucket)2003, 2007:incexpose.std:ses3"] <- c(-1,0,1)
      matrix[c(28,32,36), "factor(year_bucket)2009, 2012:incexpose.std:ses3"] <- c(-1,0,1)
      matrix[c(38,42,46), "factor(year_bucket)1992, 1997, 1999:incexpose.std:ses4"] <- c(-1,0,1)
      matrix[c(39,43,47), "factor(year_bucket)2003, 2007:incexpose.std:ses4"] <- c(-1,0,1)
      matrix[c(40,44,48), "factor(year_bucket)2009, 2012:incexpose.std:ses4"] <- c(-1,0,1)
      matrix[c(50,54,58), "factor(year_bucket)1992, 1997, 1999:incexpose.std:ses5"] <- c(-1,0,1)
      matrix[c(51,55,59), "factor(year_bucket)2003, 2007:incexpose.std:ses5"] <- c(-1,0,1)
      matrix[c(52,56,60), "factor(year_bucket)2009, 2012:incexpose.std:ses5"] <- c(-1,0,1)
      
    } else {
      matrix[13:24, "segregation.std:ses2"] <- c(rep(-1, 4),
                                                 rep(0, 4),
                                                 rep(1, 4))
      matrix[25:36, "segregation.std:ses3"] <- c(rep(-1, 4),
                                                 rep(0, 4),
                                                 rep(1, 4))
      matrix[37:48, "segregation.std:ses4"] <- c(rep(-1, 4),
                                                 rep(0, 4),
                                                 rep(1, 4))
      matrix[49:60, "segregation.std:ses5"] <- c(rep(-1, 4),
                                                 rep(0, 4),
                                                 rep(1, 4))
      ## year: segregation:ses interaction
      matrix[c(14,18,22), "factor(year_bucket)1992, 1997, 1999:segregation.std:ses2"] <- c(-1,0,1)
      matrix[c(15,19,23), "factor(year_bucket)2003, 2007:segregation.std:ses2"] <- c(-1,0,1)
      matrix[c(16,20,24), "factor(year_bucket)2009, 2012:segregation.std:ses2"] <- c(-1,0,1)
      matrix[c(26,30,34), "factor(year_bucket)1992, 1997, 1999:segregation.std:ses3"] <- c(-1,0,1)
      matrix[c(27,31,35), "factor(year_bucket)2003, 2007:segregation.std:ses3"] <- c(-1,0,1)
      matrix[c(28,32,36), "factor(year_bucket)2009, 2012:segregation.std:ses3"] <- c(-1,0,1)
      matrix[c(38,42,46), "factor(year_bucket)1992, 1997, 1999:segregation.std:ses4"] <- c(-1,0,1)
      matrix[c(39,43,47), "factor(year_bucket)2003, 2007:segregation.std:ses4"] <- c(-1,0,1)
      matrix[c(40,44,48), "factor(year_bucket)2009, 2012:segregation.std:ses4"] <- c(-1,0,1)
      matrix[c(50,54,58), "factor(year_bucket)1992, 1997, 1999:segregation.std:ses5"] <- c(-1,0,1)
      matrix[c(51,55,59), "factor(year_bucket)2003, 2007:segregation.std:ses5"] <- c(-1,0,1)
      matrix[c(52,56,60), "factor(year_bucket)2009, 2012:segregation.std:ses5"] <- c(-1,0,1)
      
      
    }
    
    
    
    
  }
  
  
  
  sim <- base::apply(beta.draws, 1, function(x){
    vector <- matrix %*% x
    return(vector)
  })
  sim <- t(sim)
  
  estimates_msas[[x]] <- estimates.fct(sim)
  
  print(x)
}

for(i in 1:length(estimates_msas)){
  estimates_msas[[i]]$year <- levels(factor(data.MSA$year_bucket))
  estimates_msas[[i]]$year <- dplyr::recode(estimates_msas[[i]]$year,
                                            `1987, 1988` = "1987",
                                            `1992, 1997, 1999` = "1995",
                                            `2003, 2007` = "2005",
                                            `2009, 2012` = "2010")
  estimates_msas[[i]]$year <- as.Date(paste0(estimates_msas[[i]]$year,
                                             "-01-01"))
  estimates_msas[[i]]$std_deviations <- c(rep(-1, 4),
                                          rep(0, 4),
                                          rep(1, 4))
  if(i %in% c(1,2)){
    estimates_msas[[i]]$group <- c(rep("Democrat", 12),
                                   rep("Republican", 12),
                                   rep("Indep/Other",12))
    estimates_msas[[i]]$comparison <- "Party"
  } else {
    estimates_msas[[i]]$group <- c(rep("ses1", 12),
                                   rep("ses2", 12),
                                   rep("ses3",12),
                                   rep("ses4", 12),
                                   rep("ses5", 12))
    estimates_msas[[i]]$comparison <- "ses"
  }
  if(i %in% c(1,4)){
    estimates_msas[[i]]$variable <- "Std. Gini"
    
  } else {
    estimates_msas[[i]]$variable <- "Std. Segregation"
  }
}

estimates_msas <- bind_rows(estimates_msas)




ggplot(estimates_msas[estimates_msas$comparison == "Party" &
                        estimates_msas$std_deviations != 0,],
       aes(x = year,
           color = factor(std_deviations),
           y = Est,
           ymin = Lower,
           ymax = Upper)) +
  facet_grid(variable ~ group) +
  geom_line() +
  geom_point(position = position_dodge(2 * 365)) +
  geom_errorbar(width = 0,
                position = position_dodge(2 * 365)) +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(x = NULL,
       y = "Estimated Percent of Respondents\nWho Think 'Rich Getting Richer'",
       color = "Standard Deviations from\nSurvey Mean") +
  scale_color_manual(values = cbbPalette[c(2,3)])

## figures/local_party_time.pdf
## 5x7

ggplot(estimates_msas[estimates_msas$comparison == "ses" &
                        estimates_msas$std_deviations != 0,],
       aes(x = year,
           color = factor(std_deviations),
           y = Est,
           ymin = Lower,
           ymax = Upper)) +
  facet_grid(variable ~ group) +
  geom_point() +
  geom_line() +
  geom_errorbar() +
  theme_bw() +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(x = NULL,
       y = "Estimated Percent of Respondents\nWho Think 'Rich Getting Richer'",
       color = "Standard Deviations from\nSurvey Mean") +
  scale_color_manual(values = cbbPalette[c(2,3)])
## figures/local_ses_time.pdf
## 5x9
